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A  Three-Dimensional  FDTD-PML  Algorithm 
Based  on  Piecewise-Linear  Approximation  for  Linear  Dispersive  Media 

S.  J.  Yakura  and  D.  Dietz 

Air  Force  Research  Laboratory,  Directed  Energy  Directorate 
Kirtland  AFB,  New  Mexico  87117 


Abstract 

Starting  with  the  unsplit-field  uniaxial  PML  formulation ,  a  second-order  accurate  FDTD-PML  algorithm  is 
obtained  using  the  piecewise-linear  approximation.  Use  of  the  FDTD-PML  algorithm  results  in  the  proper 
long  time  limit  behavior  where  the  electric  field  value  decrease  exponentially  to  zero  inside  a  PML  medium 
long  after  an  electromagnetic  pulse  is  incident  on  the  PML  medium.  The  behavior  is  consistent  with  the 
other  PML  algorithm,  such  as  Gedney's  two-step  approach. 

L  INTRODUCTION 

With  the  advent  of  high  power  computers  that  provide  fast  execution  times  and  great  quantities  of 
computer  memory,  we  are  at  the  stage  where  we  can  perform  direct  numerical  calculations  of  Maxwell’s 
equations.  Out  of  many  numerical  techniques  available  in  the  computational  electromagnetic  community, 
one  that  has  shown  a  great  promise  in  the  time  domain  is  the  well-known  finite-difference  time-domain 
(FDTD)  method  [1].  It  is  based  on  using  a  simple  staggered  differencing  scheme  in  both  time  and  space  to 
calculate  the  transient  behavior  of  electromagnetic  field  quantities.  One  of  the  greatest  challenges  of  the 
FDTD  methods  has  been  the  efficient  and  accurate  formulation  of  electromagnetic  wave  interactions  in 
unbounded  regions.  For  such  problems,  an  absorbing  boundary  condition  must  be  introduced  at  the  outer 
layer  boundary  to  simulate  the  extension  of  the  lattice  to  infinity.  One  approach  that  has  given  a  great 
promise  in  realizing  such  an  absorbing  outer  boundary  inside  the  finite  volume  computational  domain  is  the 
well-known  perfectly-matched-layer  (PML)  algorithm  that  was  first  introduced  by  J.  P.  Berenger  [2]  in 
1994  for  the  free  space  boundary  interface.  Since  that  time  Chew  and  Weedon  [3]  came  up  with  the 
modified  PML  algorithm  that  is  based  on  complex  coordinate  stretching,  which  is  shown  to  be  equivalent 
to  the  anistropic  PML  medium  approach  [4]. 

In  this  paper,  we  explore  the  formulation  of  a  3-dimensional  PML  algorithm  used  in  outer  layer 
absorbing  boundary  of  a  dispersive  medium  to  absorb  all  outgoing  waves  out  of  a  finite  simulation  volume. 
We  consider  the  case  where  a  plane  wave  propagates  outwardly  from  a  dispersive  medium  to  the  dispersive 
PML  medium  through  a  reflectionless  PML  interface.  We  start  the  analysis  based  on  unsplit-field  uniaxial 
PML  formulation  [4-8]  of  Maxwell's  equations  that  are  obtained  in  the  frequency  domain  inside  the 
dispersive  PML  medium.  We  perform  the  inverse  Fourier  transform  of  these  equations  from  the  frequency 
domain  to  the  time  domain  to  obtain  a  set  of  ordinary  first-order  differential  equations.  Then,  these 
equations  are  finite  differenced  in  both  time  and  space  using  the  usual  Yee  FDTD  scheme  while  expanding 
the  electric  and  magnetic  field  vectors  in  time  using  the  Taylor  series  expansion  about  the  current  time  step 
in  order  to  evaluate  next  time  step  values  of  the  electromagnetic  field  quantities.  Depending  on  the  number 
of  terms  kept  in  the  Taylor  series  expansion,  we  can  numerically  evaluate  the  updated  values  to  any  desired 
accuracy  we  want.  In  Section  II,  we  use  the  piecewise-linear  approximation,  which  is  equivalent  to  using 
only  the  first-order,  time-dependent  term  of  the  Taylor  series  expansion,  to  show  the  process  involved  in 
obtaining  a  second-order  accurate  FDTD-PML  algorithm.  To  obtain  higher-order  accurate  FDTD-PML 
algorithms  in  time,  we  simply  need  to  include  higher-order,  time-dependent  terms  in  the  Taylor  series 
expansion  and  follow  the  same  steps  shown  in  Section  II.  A  consequence  of  the  higher-order  accurate 
FDTD-PML  algorithm  is  the  need  to  solve  for  zeroes  of  the  nth  degree  polynomials  at  each  time  step  in 
order  to  update  field  values.  It  arises  because  of  the  use  of  the  nth-order,  time-dependent  term  of  the  Taylor 
series  expansion.  The  numerical  process  involved  in  updating  the  field  values  is  similar  to  the  FDTD 
algorithm  obtained  for  nonlinear  dispersive  media  [9,10]. 
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II.  PML  FORMULATION  FOR  LINEAR  DISPERSIVE  MEDIA 


For  a  wave  propagating  into  anisotropic,  uniaxial  dispersive  PML  media,  the  modified  Maxwell’s 
equations  under  the  PML  formulation  with  stretched  coordinates  [3]  can  be  expressed  in  the  frequency 
domain  (e'“‘  convention)  as 

V  xE(co;  x)  =  -  ico  SPML  (to)  •  p0pR  H(to;  x) , 

V  xH(00;x)  =  ico  SPML((o)*D(co;x), 
with 


(2.1) 

(2.2) 


Pfnax 

D(co;  x)  =  8o£rE(co;x)  +  £o  Pp  (co;  x)  , 

p=i 


Pp(co;x)  =  Xp!)  (co)  E(o;x), 


(23) 

(2.4) 


where  E(co;x)  is  the  electric  field  vector,  H(co;x)  is  the  magnetic  field  vector,  D(cd;x)  is  the  displacement 
field  vector,  Pp(co;x)  is  the  electric  polarization  vector,  SPML(co)  is  the  uniaxial  anisotropic  PML  matrix,  e0 

is  the  free  space  electric  permittivity,  £r  is  the  relative  permittivity,  Po  is  the  free-space  permeability,  Pr  is 
the  relative  permeability,  and  Xp(1)(co)  is  the  pth  term  of  the  collection  consisting  of  p,^  frequency- 
dependent,  first-order  (linear)  electric  susceptibility  functions,  where  p,^  is  the  maximum  number  of  terms 
which  we  choose  to  consider  for  a  particular  formulation  of  Eq.  (2.3).  Also  seen  in  the  above  equations  is 
the  notation  •  that  is  used  to  denote  a  dot  product.  Elements  of  the  uniaxial  anisotropic  PML  matrix, 
SPML  (oa),  are  given  by 

f  Sy(©)S2(co)  ) 

— -  0  0 

Sx(o» 


SPML  (co)  = 


Sx(to)Sz((o) 

Sy(CO) 


(2.5) 


0  0  Sx(m)Sy(co) 

V  j 

where  Sx(to),  Sy(co)  and  Sz(to)  are  arbitrarily  defined  co-dependent  functions  that  satisfy  the  impedance 
matching  condition  at  the  interface  of  the  non-PML  medium  and  the  PML  medium.  It  is  a  common  practice 
in  the  FDTD  community  to  choose  Sx(co),  Sy((0)  and  Sz(co)  in  the  following  forms: 


* 


Sx(co)  -  1+  °x 

with 

Ox  _ 

Ox 

(2.6-2.7) 

iCOEoER 

EoEr 

^  av 

Ov 

* 

S  (CO)  =  1  + - — 

with 

y  _ 

-  ,  and 

(2.8-2.9) 

iCOEoER 

EoEr 

fioJlR 

Sz(a»  =  l+-^_ 

with 

Oz  _ 

* 

(2.10-2.11) 

ICOEoEr 

EoEr 

where  ox  ,  ay  and  cz  are  the  PML  electric  conductivities,  and  ox ,  q*  and  Oz  arc  the  PML  magnetic 

conductivities  with  subscripts  x,  y  and  z  denoting  the  directions  in  which  PML  conductivities  are  assigned 
[2].  These  PML  conductivities  are  introduced  arbitrarily  in  order  to  implement  the  FDTD-PML  algorithm. 

We  first  eliminate  D(o);x)  in  favor  of  expressing  Maxwell’s  equations  in  terms  of  E(co;x)  and  Pp(to;x)  by 
substituting  Eq.  (2.3)  into  Eq.  (2.2).  Upon  taking  the  inverse  Fourier  transforms  of  Eqs.  (2.1),  (2.2)  and 
(2.4)  and  using  the  expressions  shown  in  Eqs.  (2.6)  through  (2.11),  we  can  show  after  some  manipulations 
Eqs.  (2.1),  (2.2)  and  (2.4)  are  written  in  the  following  time-dependent  equations: 
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(2.12) 


IH-iR  +  |a0|VPo.H(t;x)  +  Ho|XR^Fi»HDelay(t;x)  + V  xE(t;x)  =  0 , 


3t  ^  = 
I(t;  j 

at 


3E(t;x)  yi  ^Ep(l>2L)  r  Tin  .n 

BoEr - ^  ~  +  Bo^ - — - +  £q^F o  •  [£RE(t;  x)  +  2^Pp  (t;  x)] 

p  p 


+  eoTi  •  [brE06^  (t;  x)  +  J  P^lay  (t;  x)]  ~  V  x  H(t;  x)  =  0 , 


with 


oo 

Pp(t;x)=|dTX«>(t-x)E(x;x), 

^oo 

t 

jj  Delay  (t;x)  =  J*  dx  <E>(t  -  x)  •  H(x;  x) , 

— oc 

t 

EDelay(t;x)  =  Jdx£(t-x).E(x;x), 

— OO 

t 

E™ay  (t;  x)  =  Jdx  0(t -x)  *Pp  (x;  x) , 


(2.13) 

(2.14) 

(2.15) 

(2.16) 

(2.17) 


where 


^0  = 


CT 

BoBr 


-  °x 
BoBr  BoBr 


o 


°x  |  °z 
BoBr  BoBr 


BoBr 


BoBr  BoBr 


wz 

BoBr 


7; 


(2.18) 


¥1  = 


/ 

fCTy 

°z 

qx  ^ 

BoBr 

BoBr 

BoBr 

BoBr 

J 

V 

/ 

Q  ^ 

y 

f  °Z 

O  ^ 

uy 

BoBr 

V 

BoBr 

J 

BoBr 

V 

BoBr 

y 

Ox  <*z  Y  °y  °z 


BoBr  BoBr 


A 


O(t-x)  = 


BoBr 

A 


BoBr 


(2.19) 


J) 


exp[-( — — )(t — x)] 

BoBr 

0 

0 


0 

exp[-  (-— ^— )(t  -  x)] 

BoBr 


exp[-( — — )(t  —  x)] 

BoBr 


(2.20) 


In  the  above,  HDelay(t;x),  E^AUx)  and  PpDeIay(t;x)  are  introduced  to  handle  the  delayed  time-response 
behavior  of  H(t;x),  E(t;x)  and  Pp(t;x),  respectively.  These  functions  follow  naturally  from  taking  the  inverse 
Fourier  transforms  of  convolution  functions  [l/(io)i  +  A)]  H(co;x),  [l/(icoi  + A )]  E((o;x),  and 

[l/(icoi  +  A)]  Pp(co;x)  by  realizing  the  inverse  Fourier  transform  of  [l/(icoi  +  A)]  is  given  by  exp(-At), 


where  A  is  a  time  independent  diagonal  matrix  expressed  as  diag  [  ox  A^rX  ay  A^rX  az  A^oEr)]* 

To  solve  Eqs.  (2.12),  (2.13),  (2.14),  (2.15),  (2.16)  and  (2.17),  we  need  to  specify  the  expression  for  the 
linear  electric  susceptibility  function.  In  this  paper  we  consider  the  case  in  which  the  linear  electric 
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susceptibility  function  is  expressed  as  a  complex  function  that  contains  complex  constant  coefficients  and 
exhibit  exponential  behavior  in  the  time  domain  as  follows: 

X^(t)  =  Re  {apexp[-(Yp)t] }  U(t) ,  (2.21) 

where  Re{  }  is  used  to  represent  the  real  part  of  a  complex  function,  U(t)  is  the  unit  step  function,  and  ctp 
and  yp  are  complex  constant  coefficients.  Now,  Eq.  (2.14)  takes  the  form 

t 

Pp(t;x)  =Re{Qp(t;x)}  =  Re{ap  Jdtexp[-(Yp)(t-x)]E(x;x) } ,  (2.22) 

where  complex  function  Qp  (t;x)  is  introduced  in  the  above  equation  such  that  the  real  part  of  the  complex 
function  results  in  Pp(t;x). 

We  need  to  point  out  that  by  making  the  proper  choices  of  complex  constant  coefficients  and  performing 
the  Fourier  transform  of  Eq.  (2.21),  we  can  readily  obtain  the  familiar  constant  conductivity  [i.e.,apis  real 

and  yp  =  0],  Debye  [i.e.,  ap  and  yp  are  both  real]  and  Lorentz  [i.e.,  ap  is  imaginary  and  Yp  is  real]  forms  of  the 
complex  permittivity  in  the  frequency  domain. 

To  derive  FDTD  expressions  based  on  Yee  FDTD  scheme,  Eqs.  (2.12),  (2.13),  (2.14),  (2.15),  (2.16)  and 
(2.17)  have  to  be  solved  numerically  for  H(t;x),  E(t;x),  Pp(t;x),  H^^^tjx),  EDelay(t;x)  and  PpDclay(t;x)  at  each 
time  step  by  correctly  carrying  out  the  numerical  integration  of  convolution  integrals  Pp(t;x),  HDeIay(t;x), 
EDe'ay(t;x)  and  PpDelay(t;x).  Therefore,  the  whole  solution  rests  on  the  question  of  how  to  carry  out  the 
numerical  integration  of  Pp(t;x),  HDclay(t;x),  EDelay(t;x)  and  PpDelay(t;x)  at  each  successive  time  step.  For  that 
reason,  the  rest  of  this  section  is  devoted  to  the  numerical  formulation  that  treats  Pp(t;x),  HDelay(t;x), 
EDeiay(t;x)  and  Pp^^tyc)  into  the  overall  FDTD  scheme  based  on  the  recursive  convolution  approach. 

We  first  convert  the  convolution  integrals  Pp(t;x),  HDclay(t;x),  EDelay(t;x)  and  PpDeIay(t;x)  [i.e.,  Eqs.  (2.22), 
(2.15),  (2.16)  and  (2.17)]  into  the  following  equivalent  first-order  differential  equations: 


9Qp(t;x) 


gt  +(yP)Qp(t;i)  =  apE(t;x) , 

(2.23) 

agPegyt(t;-)  +£(t).HDel^(t;x)  =  H(t;x), 

(2.24) 

~g  +£(0*EDelay(t;x)  =  E(t;x), 

(2.25) 

3Q°elay(t;x) 

- +£(0*Q^  y(t;x)=Qp(t;x), 

(2.26) 

where  complex  function  Qp1*133'  (t;x)  is  introduced  in  Eq.  (2.26)  such  that  the  real  part  of  the  complex 
function  results  in  PpDelay(t;x). 

To  show  how  we  can  use  Eqs.  (2.12),  (2.13),  (2.23),  (2.24),  (2.25)  and  (2.26)  to  come  up  with  a  3-D 
FDTD-PML  algorithm  for  dispersive  PML  media,  we  integrate  Eqs.  (2.12)  and  (2.24)  from  t=(n-V£)At  to 
t=(n+Vi)At,  and  Eqs.  (2.13),  (2.23),  (2.25)  and  (2.26)  from  t=nAt  to  t=(n+l)At.  Then  Eqs.  (2.23),  (2.24), 
(2.25)  and  (2.26)  are  solved  exactly  using  the  integrating  factor  technique  for  a  given  discrete  time  interval 
to  go  forward  in  time  by  At.  The  result  is  that  we  need  to  evaluate  definite  integrals  appearing  in  the 
following  equations: 

(n+'^)At  'jTT/  .  (n+fc)At 

J  dx  ~"~d  ~ ~  +  JiopR^Po.  J  dx  H(x;  x) 

(n-'/i)At  (n-lA)At 

(n+*/4)At  (n+'/i)At 

+  JdxH^Cxjx)  +  Jdx  V  xE(x;x)  =  0,  (2.27) 

(n-'^)At  (n-Vi)At 
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(n+1)Al  ac/  x  0>+DAt  30  (x‘  x) 

“>  1  j  Ot-^- 

nAt  P  nAt 


Vnri;zai  tn+ijZM 

+  eoeR^o  •  JdxE(x;x)  +eovFo«Re  JdxQp(x;x) } 

nAt  P  nAt 

(n+l)At  (n+l)At 

+  eoeR^i  •  J dx  (x;  x)  +  Eo¥i  .  Re{£  Jdx  Q™ay  (x; x) } 

nAt  P  nAt 

(n+l)At 

-  Jdx  VxH(x;  jc)  =  0,  (2 

nAt 

(n+l)At 

Qp  (nAt + At;  x)  =  exp[-(Yp)At]  [  Qp  (nAt;  x)  +  ctP  J  dx  exp[-(Yp)(nAt  -  x)]  E(x;  x)],  (2 

nAt 

HDelay  (nAt  +  Vi  At;  x)  =  exp(-OAt)  •  [  HDelay  (nAt  -  ViAt;  x) 

(n+'/4)Al 

+  J  dx  exp[-4>(nAt  -  Vi  At  -  x)]  •  H(x;  x)  ] ,  (2 

(n-H)At 

(n+l)At 

E  Delay  (nAt + At;  x)  =  exp(-4>At)  •  [  E0*1^  (nAt;  jc)  +  J  dx  exp[-4>(nAt  -  x)]  •  E(x;  x)  ]  ,  (2 

nAt 

(n+l)At 

q Delay  +  ^  _  eXp(-OAt)  •  [  (nAt;  x)  +  J  dx  exp[-0(nAt  -  x)]  •  Qp  (x;  x)  ] 

nAt 

(n+l)At 

=  exp(-OAt)  •  [  Q^lay  (nAt;  x)  +  J  dx  exp[-£(nAt  -  x)  - 1  (Yp)(x  -  nAt)]  •  Qp  (nAt;  x) 


(n+l)At  t 


+  aP  J  dx  J  dx’  exp[-<t>( nAt  -  x)  - 1  (Yp)(x  -  x’)]  •  E(x’;  x)  ]  ,  (2.32) 

nAt  nAt 

where  I  is  the  identity  matrix.  Furthermore,  some  of  those  definite  integrals  that  appear  in  Eqs.  (2.27)  and 
(2.28)  are  manipulated  and  cast  in  the  following  forms: 

(n+‘/2)At  (n+»/2)At 

J  dx  H0*1^  (x;x)  =  Jdx  exp[-4>  (x-  nAt +14  At)]  •  HDelay  (nAt-V4At;x) 

(n-^)At  (n->/2)At 

(n+‘>4)At  x 

+  Jdx  Jdx’exp[-|>(x-x’)]*H(x’;x) ,  (2.33) 

(n-'/2)At  (n-V£)At 


\UT 

J  dxQp(x;x)  =  J  dxexp[-(Yp)(x-nAt)]Qp(nAt;x) 


(n+l)At  x 


■  ap  J  dx  J  dx’exp[-(Yp)(x  -  x’)]  E(x’;x) , 
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(n+l)At 


(n+l)At 


JdxE^'^^x)  =  j  dx  exp[-4>  (t  -  nAt)]  •  E  Delay  (nAt;  x) 

+l)At  T 

JdT  JdT’exp[-<b(T-T’)]»E(T’;x)> 

nAt  nAt 

J dT  Q^lay  (t;  x)  =  JdT  exp[-d>  (t  -  nAt)]  •  Q^lay  (nAt;  x) 

nAt 

(n+l)At  x 

+  JdT  JdT’exp[-£T-I(Yp)(T-nAt)]*  Qp(nAt;x) 


nAt 


(n+l)At 


(2.35) 


(n+l)At 


nAt  nAt 
(n+l)At 


+  aP  JdT  JdT’JdT"exp[-£T-I(Yp)(T’-T")]*E(T";x).  (2.36) 

nAt  nAt  nAt 

To  obtain  second-order  accuracy  in  time  from  a  finite  differencing  technique,  H(t;x)  and  E(t;x)  are  taken 
to  be  piecewise-linear  continuous  functions  over  the  entire  temporal  integration  range  such  that  H(t;x)  and 
E(t;x)  change  linearly  with  respect  to  time  over  given  discrete  time  step  intervals.  It  is  equivalent  to  using 
only  the  first-order,  time-dependent  terms  of  the  Taylor  series  expansions  for  H(t;x)  and  E(t;x), 
respectively,  expanded  in  time  about  current  time  step  of  H(t;x)  and  E(t;x).  Mathematically,  we  can  express 

H(t;x)  and  E(t;x)  in  the  following  forms  in  terms  of  CtDiji."  ^,  (tI)iji,n+!/;,  (E)ijkn  and  (E)ijkn+I  where  superscripts 
n-'/2,  n,  n+'A  and  n+1  are  used  to  denote  discrete  time  steps  at  t=(n-!4)At,  t=nAt,  t=(n+lA)At  and  t=(n+l)At, 
respectively.  Subscripts  are  used  to  denote  discrete  spatial  locations,  x=[iAx,  jAy,  kAz]  for  E(t;x)  and 

x=[(M£)Ax,  (j-V4)Ay,  (k-Vi)Az]  for  H(t;x),  with  Ax,  Ay  and  Az  being  the  spatial  grid  sizes  in  the  x,  y  and  z 
directions,  respectively. 


H(l;x)  = 


At 


0, 


-[t-(n  -Vi)At]  + higher  order  terms, 

for  0  <  (n  -  Vi)  At  <  t  <  (n  +  Vi)At  (2.37) 

for  t  <  0 


l(t;x) 


_n  ,  [(E)£*-(E)Sk] . 

(E)^  + - — - (t  -  nAt)  +  higher  order  terms, 


for  0  <  nAt  <  t  <  (n  +  l)At 


(2.38) 


(0>  for  t  <  0 

By  keeping  more  than  the  first-order,  time-dependent  term  in  the  above  Taylor  series  expansion,  it  is 
possible  to  investigate  higher  than  second-order  accurate  FDTD  algorithms. 

Substituting  Eqs.  (2.37)  and  (2.38)  into  Eqs.  (2.27)  through  (2.36),  performing  the  time  integration  from 
t=(n-!4)At  to  t=(n+Vi)At  for  field  values  that  depend  on  the  magnetic  field  [i.e.,  H(t;x)  and  HDelaj'(t;x)],  and 
from  t=nAt  to  t=(n+l)At  for  field  values  that  depend  on  the  electric  field  [i.e.,  E(t;x),  Qp(t;x),  EDelay(t;x)  and 
Qp  y(t;x)]  and  using  the  usual  staggered  Yee  FDTD  scheme  for  spatial  discretization,  we  obtain  a 
3-dimensional  FDTD-PML  algorithm  in  Cartesian  coordinates.  .  After  some  manipulations,  we  obtain  the 
following  six  equations  which  constitute  the  entire  set  of  updating  expressions  needed  to  update  field 

values  (E)y1(n+1,  (H)ijkn+I,  (Qp)ijt"+\  (EDelay)ijkn+1,  (HDe,ay)ijk”+1  and  (QpDelay)ijt“+1  inside  a  dispersive  PML 
medium: 


Qo  .(H)Hk+'/j  +Q,.(H)^  +£i2.(HDe,ay)H^  +Se  =0, 


(2.39) 
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+2,.®"*%  +  £i.Re{X(Q“”);)+s!1  =  0, 

P 

(Qp)ijk1  =0P,o  [(Qp)."k  +©P,i  (Dijk  +©P,2  (Dp1], 

(HDelay)lT  =  £6  •[  (HDelay)1^  +07  •(H)Jj^  +£8*(H)p/2]. 

(EDelay)^1  =£6«[(EDe,ay)(;k  +£7 -(Dp  +£g*(D.“k1], 

(QpDelay)fjk+1  -£6-[(QpDelay)5k  +n^o*(Qp)1"k  +£p,  .(E)sk  +np2«(E)^1], 

where  SE  and  SH  are  given  by 

At  r  ait-  v  n  A-r-  v  n  T  At 


Sp  “ 


(^)Ay[(Ez),“(i+,'4)k  (Ez)'Wi  (Wr)Az 


At 


[(E  x  )  jj(k+V4)  (Ex  ] 


At 


(Ho|Jk)Az 

(HoM*)Ax[(Ey)“i^)jk  ” (Ey)(i-*)jk  J  (^)Ay 


[(Ey^-CEy),^] 

[(E  z )  (!+t/2)jk  “  (E  z )  (S->/2)ik  ] 


ij(k^)J  (^)AxLV^^(i+1^ 

ik  ] “ ~  z  x  a  [(^ X  )f(j+,/2)k -(Ex )“a^)k ] 


z  /  (i-*/2)jk . 
n 

i(j — V^)k  - 


sH  = 


At  [(Hz ~ (Hz)^)k ] +_^L-[(Hy )1^+,/!) -(H,)^] 


-(Hz) 


n+»/2 


z  /  (iH6)jk 


(&>6i)Ay~  (E.oiaz 

»•  -<H» 


(2.40) 

(2.41) 

(2.42) 

(2.43) 

(2.44) 


(2.45) 


(2.46) 


(EoeR)Ax 

\ 

and  0pO,  ©pl  and  ©p2  are  the  time-invariant  coefficients,  and  Q0  ,  ,  Q2  ,  Q3  ,  £24  ,  Q5  ,  Q6  , 

Q7  ,  Q.%  ,rp0  ,  np0  ,  npl  and  np2  are  the  time-invariant  matrices.  Both  time-invariant  coefficients 
and  matrices  depend  only  on  material  properties  aP ,  Yp ,  ax ,  ay  and  oz »  and  time  increment  At.  Shown 


in  Appendix  A  are  the  explicit  expressions  of  these  coefficients  and  matrices  expressed  in  terms  of  aP  *  Yp » 
Ox »  Oy »  az  and  At. 

Using  the  above  FDTD-PML  algorithm  the  computer  simulation  can  be  performed  for  electromagnetic 
waves  that  propagate  inside  dispersive  PML  media  by  simply  going  through  the  following  steps: 

(1)  First,  before  updating  the  field  values,  time-invariant  coefficients  0p  O,  ©p  l  and  0p  2 ,  and 

time-invariant  matrices  fl0  ,  Qx ,  Q2  ,  Q3  ,  Q4  ,  Qs  ,  fi6  ,  Q7  ,  fl8  ,rp0  ,  np0  ,  np  l 
and  np  2  are  calculated  at  the  beginning  of  the  simulation  as  part  of  the  initial  condition  for 


given  values  of  ap  >  YP  >  ax  *  ay  *  az  and  At.  These  calculated  values  are  stored  in  memory 
and  used  at  each  time  step  to  update  field  values  (E)ijkn+1,  (H)ijkn+1,  (Qp)ijkn+1,  (EDeIay)ijkn+l, 
(HDelay)ijkn+1  andfQ^'V. 

(2)  Using  Eq.  (2.39),  (H)ykn4%  is  calculated  based  on  the  known  values  of  (HW1’^  and  (HDelay)ijkD  !4 
and  (E)ijk“ . 

(3)  Using  Eq.  (2.42),  (HIX'lay)jjkn+,/;  is  calculated  based  on  the  known  values  of  (HDelay)ykII‘ys, 
(H)ijkn+/2  and  (H)ykn 

(4)  Using  Eq.  (2.40),  (E)ykn+1  is  calculated  based  on  the  known  values  of  (E)ykn,  (EDelay)yk11, 
(QP)ijkn,  (Qp^ijk"  and  (H)ijkn+“  . 


(5)  Using  Eqs.  (2.41),  (2.43)  and  (2.44),  (Qp)ijkn+1,  (EDe,ay)ijln+l  and  (Qp^V*'  are  calculated 
based  on  the  known  values  of  (E)ijk"+1,  (E)ijk°,  (Qp)ijkn,  (EDe,ay)ijkn  and  (QpDelay)ijkn . 

(6)  Increment  the  time  step  by  At.  Go  back  to  step  (2)  and  repeat  the  whole  process  over  again. 
Shown  in  Figure  1  is  the  flow  chart  of  numerical  steps  required  to  update  field  values  as  described  above. 

In  the  case  of  the  PML  interface  to  the  constant  conductivity  medium  we  setyp  =  0  and  ap  to  be  the 
value  of  the  constant  electric  conductivity.  The  result  is  that  the  matrix  elements  simplify  to  the  forms 
shown  in  Appendix  B.  Furthermore,  if  we  set  both  apand  yp  to  zeroes  the  FDTD-PML  algorithm  reduces 
to  the  case  of  the  simple  PML  algorithm  for  the  vacuum. 

III.  CONCLUSIONS 

We  present  in  this  paper  the  formulation  of  a  three-dimensional  FDTD-PML  algorithm  inside  dispersive 
PML  media  that  is  used  to  absorb  all  outgoing  electromagnetic  waves  within  a  finite  simulation  volume  to 
create  the  notion  of  infinity  at  the  outer  layer  boundary  of  the  computational  volume.  Because  of  the  use  of 
the  piecewise-linear  approximation,  the  FDTD-PML  algorithm  provides  second-order  accuracy  in  time  for 
the  calculation  of  electromagnetic  field  quantities  inside  the  outer  absorbing  layer  boundary,  resulting  in 
less  than  one  thousandths  of  the  outgoing  wave  coming  back  into  the  main  computational  volume. 
Computationally,  we  can  see  that  the  FDTD-PML  algorithm  retains  all  the  advantages  of  the  usual  first- 
order  discrete  recursive  convolution  approach,  such  as  fast  computational  speed  and  efficient  use  of 
computer  memory.  In  the  limit  of  no  PML  interface,  the  FDTD-PML  algorithm  reduces  to  a  simple  FDTD 
algorithm  formulated  for  linear  dispersive  media. 

Lastly,  we  need  to  point  out  that  the  exponential  form  of  the  susceptibility  function  is  crucial  in  allowing 
us  to  implement  the  recursive  feature  in  our  algorithm. 
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APPENDIX  A 


This  appendix  gives  the  explicit  expressions  of  coefficients  and  matrices  seen  in  Eqs.  (2.39)  through 
(2.44).  The  coefficients  are  ©p0,  ©Pi]  and  @p  2  •  The  matrices  are  £20  ,  Q,  ,  Q2  »  G3  >  G4  »  Gs  .  G6  » 

£27  ,  Qg  ,  Tp  q  ,  np0  ,  np ,  and  np  2  .  Also,  to  express  these  coefficients  and  matrices  in  more  compact 

forms,  additional  terms,  such  as  £p0,  £pl,  £p2,  <;pl,  <;p2,  (t|/0)x,  (\p,)x ,  (v|/2)x,  ((p,)x,  (<p2)x, 

(Cp>0)x  >  (Cpj)x  >  (jtp.o)  x  >  (Jtp,i)x  an<^  (itp,2)x  ’  are  defined.  These  additional  terms  are  shown  following  the 
expressions  for  coefficients  and  matrices. 


®p,o  £p,o « 

e»..  -^itw-twl. 

8r 

«  -ap  p 

Op,2-— ^2. 


Qn  — 


(Go) 


V 


oni 

0 

0 


0 

(Go)  22 
0 


0 

0 

(gq)33 


with  three  diagonal  elements  expressed  as 


with  three  diagonal  elements  expressed  as 


(Go )n  =  1+ [(— )  +  (— )  -  (— )]—+ K— )  -  (—)][(—)  -  (—)]  (cp2) x . 

EoEr  £o£r  EoEr  2  EoEr  £o£r  EoEr  EoEr 

(Go) 22  =  Replace  [x— »y,  y-»z,  and  z— »x]  in  (Q0)n , 

(Q0)33  =  Replace  [x-»z,  y-»x,  and  z-»y]  in  (Q0)u , 

(Gj)ii  0  0 

G^  =  0  (Q,)22  0 

I  0  0  (G,  )33 

V  / 

(G, ) , ,  =  -I + [(-^)  +  -  (-^-)]  ^ + [(-^-)  -  (-^)]  [(-^-)  -  (-^)]  [(<p,), 

EoEr  EoEr  EoEr  2  EoEr  EoEr  EoEr  EoEr 
(Qt)22  =  Replace  [x— >y,  y— >z,  and  z-»x]  in  (Si^n , 

(Gi  )33  =  Replace  [x— >z,  y— »x,  and  z-»y]  in  (£2,)u , 

^(G2)„  0  0 


Qo  = 


0 

0 


(G2) 


2  222 


0 


o  (G2)33 


with  three  diagonal  elements  expressed  as 


(G2)n  =[(^L)-(^-)][(^)-(^)](Vi)x, 

EoEr  EoEr  EoEr  EoEr 

(Q2)22  =  Replace  [x-»y,  y— >z,  and  z— >x]  in  (Q2)n , 
(G2)33  =  Replace  [x— »z,  y— >x,  and  z— »y]  in  (Q2)n , 


(A.1) 

(A.2) 

(A.3) 

(A.4) 

(A.5) 

(A.6) 

(A.7) 

(A.8) 


-(tp2)x],  (A.9) 

(A.  10) 
(A.11) 


(A.  12) 


(A.  13) 


Q?  = 


(G3)n 

0 

0 


0 

(G3)22 

0 


0 
0 

(G3)33 


with  three  diagonal  elements  expressed  as 


(A.  14) 
(A.  15) 

(A.  16) 
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(A.  17) 


(Q3)n=l+[A)  +  (^)-(^)}]^+[(-^)-(-^-)][(-^-)-(-^)]((p2)> 

8o£r  £o£r  £o£r  2  £o£r  £o£r  £o£r  £o£r 


p  p 

+  7-  [(— )  -  (-^-)]  [(-^)  -  0^-)]  Re  { V  ap  (C  2)  x  } , 

£r  £o£r  £o£r  £o£r  £o£r  "  P  p,Z 


£^4  — 


(£23  )22  =  Replace  [x~>y,  y-»z,  and  z->x]  in  (Q3)n  , 
(Q3)33  =  Replace  [x->z,  y->x,  and  z-»y]  in  (Q3)u  , 
^4)11  0  0 

0  (£2  4 )  22  0 

0  0  (£24)33 

a 


with  three  diagonal  elements  expressed  as 

)][(<Pi)x  -(CP2)> 


(£24  )n  =  -1 +K-^)  +  (-^-)-(-^)]^+[A-)  -  (-^-)][(-^-)  - 

£o£r  £o£r  £o£r  2  £o£r  £o£r  £o£r  £o£r 


P  p 

+ & [  t }  -  ( e^r)]  [( e^r)  ~ ( eir)]  Re  S  “p  [(w> *  "(y«^ 


(^4)22  =  Replace  [x— >y,  y— »z,  and  z-»x]  in  (£24)n , 
(£24)33  =  Replace  [x->z,  y-*x,  and  z->y]  in  (£24)n , 


Qs  — 


(^2)11 

£r 

0 
0 


0 


£r 


-(Q2) 


'2-' 22 


0 


£r 


0 

0 

■(Q2) 


33 


q6  — 


q7  = 


£2o  — 


(Vo)*  0  0 

0  (\jr0)y  0 

0  0  (V0)  2 

/ 

[(Vi)x-(V2)J  0 

0  [(Vl)y-(V2)y] 


0 

0 


0 


0 


[(Vi)z-(V2)J 


(V2) 


2'x 


0 


0 


O  (V2)y  0 

0  0  (V|/2)z 


rp>0 


(rp,o)n  0  0 

0  (Cp,0 )  22  0 

0  (^p,0  )  33 


0 


with  three  diagonal  elements  expressed  as 


(A.  18) 
(A.  19) 

(A.20) 

I 

-9p,21  1 

(A.21) 

(A.22) 

(A.23) 

(A.24) 

(A.25) 

(A.26) 

(A.27) 

(A.28) 
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(rp.o)n  =Up.0  -11+1  (—)  +  (—)-(—)  ]  Sp, 

H  EoEr  £o£r  £o£r 

+  [(— )  -  (—)][(—)  -  (— )]  (Cp  o)  X  . 
£o£r  £o£r  £o£r  £o£r  k’ 


(rp,o>22  =  Replace  [x-»y,  y-»z,  and  z-»x]  in  (rp0)n , 
(rp,o)33  =  Replace  [x-»z,  y->x,  and  z->y]  in  (rp0)u , 


np,0  — 


“p,  -V 

~lnp,olx 

CR 

0 

0 


OCp 

Er 


(np,o)j 


0 

0 


o  —  (np,o)2 

£r 


^p,l  — 


ap 

£r 


f  (jtp,i)  x  (ftp,2)  x  1 


0 

0 


aP 

£r 


[(jlp,l)y  “(7lP,2)y  1 


0 

0 


aP 

Er 


[(7Cp<i)  z  (7Cp,2)  z 


np,2  “ 


aP  ,  , 

~1  v7Cp,2/  x 

Er 


0  —  (np,2)y 

Er  1 


0 

0 


0 


0  —  (np,2>2 

Er 


and  matrix  elements  are  defined  as  follows 
^P,o  s  exp[-(Y  p ) A  t] , 

rl— exp[-(Yp)At] 


At  , 

£P,i  s  JdT  exp[-(Yp)T]  =  At  [- 


0 

At 


SP,2  s  T dx  (77 )  exp[-(YP)T|  =  At — ^— [ 

i  At  p  <v  Ut 


(Yp)At  ^  ’ 

1  r  1  — exp[— (Yp)At] 


0 

At  T 


(Yp)At  (Yp)At 


-exp[-<YD)At]] . 


l-exp[-(YD)At] 


SP,i  =  fdx  f  dx’  exp[-(Yp)(x  -  x’)]  =  (At)2  — - — [l - —  wp/ 

J  J  P  (YP)AtL  (Yp)At 


]. 


0  0 
At  t 

Sp,2  =  JdT  J  dT’  )  exp[-(Yp)(x-x’)] 
0  0 

=  (At)2— !— [i - i — [l- 

(Yp)At  L2  (Yp)At 

For  the  x  component: 

(v0)x  sexp[-(— )At], 

EoEr 


l-exp[-(Yp)At] 

(Yp)At 


]]. 


(A.29) 

(A.30) 

(A.31) 

(A.32) 


(A.33) 


(A.34) 


(A.35) 

(A.36) 

(A.37) 

(A.38) 


(A.39) 

(A.40) 
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(A.41) 


At 

(V,)x  =  fdx  exp[-(— )x]  =  At  [ 

J  £o£r 


1 — exp[-( — — )At] 


£o£r 


(— )At 

£o£r 


At 


7xo  1  1 — exp[-( — — )At)] 

(V2)x  =  J  dT (— )  exp[-( — — )x]  =  At - [ - — - expH-^-)At)]] , 

~  At  EoEr  ,  (Tv  . .  gy 

0  (— )At  (— ^-)At 

£o£r  £o£r 

a. 


£o£r 


At  x 


(<Pi)x  =  fdx  f  dx’exp[-(-^-)(x-x’)]  =  (At)2 - - - [l- 

i  J  EoEr  ov  . 


l-exp[-( — —  )At] 


£o£r 


0  0 

At  T 


(— )At 

Eo£r 


(— )At 

£o£r 


(<P2) x  s  1 dt  [ dT’  (7- )  exp[-(— )(T  - X’)] 
i  i  At  EoEr 


0  0 


= <it)2  — — ' — [i — - — [1 

<-^)A,  2  (i-)A, 

Eo£r  £o£r 


l-expH^-)At] 


EoEr 


(— )At 

EoEr 


-]]. 


At  x 


(Cp,0)x  =  JdT  j  dx’exp[-(Yp)T’]  exp[-(^-)(x-x’)] 


0  0 


=  (At)2 


1 _ [~l-exp[-(Yp)At]  i-exp[-(^-)At] 


K-^)At-(Yp)At] 

EoEr  v 


(Yp)At 


(— )At 

EoEr 


]. 


At  x  x’ 


(Cp,i)x  =  fdx  J dT> j* dT"  exP[-(Yp)('t“T" )] exp[-(-^-)(x - x’)] 

i  \  *  EoEr 


0  0 


=(At)J — !— {  [ — ! — ][, 


l-exp[-( — — )  At] 


EoEr 


(Yp)At  \_°x_ 


4 


(— )At  (— -2!-)At 

EoEr  EoEr 


-[ 


1 


(— )At-(Yp)At 

EoEr  H 


][ 


l“cxP[~(Yp)All 


(YD)At 


(— )At 

EoEr 


■]}. 


At  x  x’ 


tCp,2)x  =  |dT  J  dT’  [  dx"  (-£- )  exp[-(y  )(x-x")]  exp[-(-^-)(x-  x’)] 

i  i  i  At  M  EoEr 


0  0 


=  (At)3 


1  rr  1  in  i-exp[-(— *-)At] 

— ]  [}-r> — — ]C^ — +—!—]] 

EoEr 


(Yp)At 


(— )At 

EoEr 


(— )At 

EoEr 


J  1  1  [  1  1  rl-exp[-(Yp)At]  1  exPt  (eoeVAt^-|  1 

L  a„  J  L/V  \A.JL  o,  n  J  /’ 


(— )At-(Yp)At  UP 

EoEr  M 


(Yo)At 


(Yp)At  ,<*x 


( — ~)At 

EoEr 


(A.42) 


(A.43) 


(A.44) 


(A.45) 


(A.46) 


(A.47) 
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f 


* 


At 

(7tp,o)x  s  Jdx  exp[-(Yp)x]exp[(— )t] 


=  Atexp[(-^)At][ 


exp[-(Yp )  A  t]  -  exp[-( — — )  At] 


EoEr 


EoEr 


(A.48) 


( — — )At-(Yp)At 


£o£r 


At  T 

(7tp.i)x  =  fdx  fdT’exp[-(Yp)(x-T’)]exp[(— )t] 

J  J  p  £o£r 


=  (At)' 


exp[( — — )At]  f  _  1  -  exp[— ( — — )At] ..  f  exp[-(Yp  )A  t]  -  exp[-(^-)At]  -  . 
EoEr  J  I  EoEr  |  |  _ _ _ EoEr _ J  I 


(Yp)At 


:{[■ 


£o£r 


(— )At 

£o£r 


( — — )At-(Yp)At 

c„c„  r 


£o£r 


At  T 


(np,2>  x  s  fdx[dT’(-^-)  exp[-  (yp)(x  -  t’)1  exp[(— )t] 

*  J  At  K  EoEr 


(A.49) 


=  (At)2 


exp[( — — )At] 

_ £o£r _ 

t(YP)At][(— 5-)At] 


M 


l-exp[-<-i-)At] 


£o£r 


£o£r 


gv 

(— )At 

£o£r 


][ 


1+ 


(— )A- 

EoEr 

(Yp)At' 


-[ 


exp[-(Yp)A  t]-  exp[-(— )At]  (— )  At . 


£o£r 


]W}, 


(Yp)At 


(A.50) 


( — — )At  “(Yp)At 

EoEr  y 

For  the  v  component: 

Replace  [x — >y]  of  matrix  elements  defined  for  the  x  component  above. 

For  the  z  component: 

Replace  [x— >z]  of  matrix  elements  defined  for  the  x  component  above. 

As  noted  in  the  definition  of  matrix  elements  above,  these  elements  depend  only  on  known  values  At,  aP  * 
Yp .  Ox  y  a y  and  az  * 


APPENDIX  B 


By  letting  yp  —>  0  in  Eqs.  (A.35)  through  (A.50)  we  can  obtain  the  following  matrix  elements  for  the  case 
of  the  constant  conductivity. 


^p,0  “  1  * 

(B.l) 

^P,i  =  At , 

(B.2) 

-  At 

Sp,2-y> 

(B.3) 

• 

.  (At)2 

Spj  2  ’ 

(B.4) 

r  (At)2 

^P,2  6 

(B.5) 
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For  the  x  component: 

(Vo)*  =expH~)At], 

Eo£r 


(Vi), 


-m[ 


1-expH— — )At] 


EoEr 


ov 

(— )At 

EoEr 


]. 


(V2)x=At— - - [l- 


1-expH— *-)At)] 


EoEr 


(— )At 

EoEr 


(— )At 

EoEr 


] 


r  l-exp[-( — — )At] 

«p,)x  =  (At)2 - - - [l - §2§£ - 1 

fT  L  r r  J 


(— )At 

EoEr 


Oy 

(— )At 

EoEr 


(<p2)x  =  (At)2  —1 - [I — _i - [l- 


1— exp[-( — — )At] 


EoEr 


(^)At  2  (^)At 

EoEr  £o£r 


(— )At 

EoEr 


-]]. 


(Cp,o)x  =(At)2— — i - [l- 


1  ~  exp[-( — — )At] 


EoEr 


-1. 


(— )At  (— *-)At 

EoEr  EoEr 


[l- 

( — — )At  2  (-^-)At 

EoEr  EoEr 


l-exp[-(-^-)At] 


EoEr 


(■~A— )At 

EoEr 


-]]. 


<w>=(40>— 


1-expH— M  At] 


EoEr 


(— *-)At 

EoEr 


EoEr 


aY 

(— )At 

EoEr 


■]]}■ 


(%>,o)x  =  At  [ 


1— exp[-( — — )At] 


EoEr 


(— )At 

EoEr 


]. 


(Jtp,i)x  =  (At)2 - - - [l- 


l-exp[-( — — )At] 


EoEr 


(-— *~)At 

EoEr 


aY 

(— )At 

EoEr 


1  —  exp[— (-^-)At] 


EoEr 


a¥ 

(— )At 

EoEr 


}l 


<*P.2>*  =  (At)2  — - ^ - [  1- 

(•^MAt  2  (— *-)At 

EoEr  £o6r 

For  the  v  component: 

Replace  [x — »y]  of  matrix  elements  defined  for  the  x  component  above, 
For  the  z  component: 

Replace  [x— >z]  of  matrix  elements  defined  for  the  x  component  above. 


(B.6) 

(B.7) 


(B.8) 


(B.9) 


(B.IO) 


(B.ll) 


(B.12) 


(B.13) 


(B.14) 


(B.15) 


(B.16) 
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Flow  Chart 
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